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1 Abstract 

We show that Bose-Einstein condensates in a honeycomb optical lattice are described by a nonlinear Dirac equation in the long 
wavelength, mean field limit. Unlike nonlinear Dirac equations posited by particle theorists, which are designed to preserve the 
principle of relativity, i.e., Poincare covariance, the nonlinear Dirac equation for Bose-Einstein condensates breaks this symmetry. 
We present a rigorous derivation of the nonlinear Dirac equation from first principles. We provide a thorough discussion of all 
symmetries broken and maintained. 
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1. Introduction 

Recently the first truly two-dimensional (2D) solid state 
material, graphene, was created in the laboratory |H2j . One 
of the most exciting aspects of this novel material is that 
long wavelength excitations are described by a Dirac equa- 
tion for massless particles, with a "speed of light" equal to 
the Fermi velocity vp ~ c/300 [3]. Thus one can study rela- 
tivistic phenomena at very low velocities in an experiment 
far more accessible than a particle accelerator. The only real 
requirement to obtain this equation is the simple hexago- 
nal, or honeycomb lattice structure of the graphene [4]. One 
can therefore consider any solid state system constructed 
on a honeycomb lattice, including artificial systems, in or- 
der to study relativistic phenomena in novel materials ac- 
cessible in tabletop experiments. 

The most precise, cleanest, most controllable artificial 
solid state system is ultra-cold atoms in optical lattices. 
Such systems have no impurities and no disorder unless 
specifically added in by hand. They are very versatile: 
they can be constructed with an arbitrary lattice structure 
in one, two or three dimensions; they can contain bosons 
and/or fermions, atoms and/or diatomic molecules; and 
they can even have a pseudo-spin structure. Their tem- 
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perature and interaction sign, magnitude, and symmetry 
can be controlled externally. Moreover, 2D physics has 
recently been of great interest in this context, in the form 
of the Berzinskii-Kosterlitz-Thouless crossover [5] , and 2D 
systems underpinned by lattices are immediately available 
in experiments. Instead of considering ultra-cold fermions, 
which could be used to produce an near-exact analog of 
graphene [BJ, we consider ultra-cold bosons. Bosonic statis- 
tics and interactions lead to a new feature in the massless 
Dirac equation known in graphene - a naturally occurring 
nonlinear term, giving rise to a nonlinear Dirac equation. 

The study of nonlinear phenomena in ultra-cold atoms, 
especially in Bose-Einstein condensates (BECs) |7I8] , 
has been enormously fruitful. The recent text edited by 
Kevrekidis, Frantzeskakis, and Carretero-Gonzalez pro- 
vides an excellent summary of this field [9] . The nonlinear 
mean field description given by the nonlinear Schrodinger 
equation (NLSE) has been very accurate in the majority 
of experiments on BECs. Vector and non-local generaliza- 
tions of the NLSE have also proven useful [9]. In optical 
lattices, the mean field description remains accurate pro- 
vided the lasers creating the standing wave which is the 
optical lattice are not too intense, and the dimensionality 
of the system is greater than one |10lllj . 

In this article, we present a completely new class of non- 
linear phenomena in BECs, based on the nonlinear Dirac 
equation (NLDE). Nonlinear Dirac equations have a long 
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history in the literature, particularly in the context of par- 
ticle and nuclear theory [12ll3ll4fL5] . but also in applied 
mathematics and nonlinear dynamics [16117118119120] . As 
nonlinearity is a ubiquitous aspect of Nature, it is natural 
to ask how nonlinearity might appear in a relativistic set- 
ting. However, this line of questioning has been strongly 
constrained by modeling, rather than first principles. That 
is, there is no standard first principle of quantum electro- 
dynamics (QED) which is nonlinear. So, the approach has 
been to require symmetry constraints in nonlinear mod- 
els. One of these constraints is the principle of relativity, 

1. e., Poincare covariance. Poincare symmetry includes ro- 
tations, translations, and Lorentz boosts. 

In contrast, our NLDE is not a model: it is derived from 
first principles for a weakly interacting bosonic gas in the 
presence of a honeycomb optical lattice. We show that 
Poincare symmetry is naturally broken by the nonlinearity 
inherent in this system. Given that this form of nonlinear- 
ity, which depends only on the local condensate density, is 
one of the most common throughout nature, it is important 
to recognize that the principle of relativity may be broken 
by small nonlinearities even at a fundamental level, for ex- 
ample of QED |21 22 . Thus, we suggest a new direction 
of investigation in particle physics of possible nonlineari- 
ties as well as providing a natural context in artificial solid 
state systems for the NLDE. 

This article is outlined as follows. In Sec. [2] we provide 
a rigorous derivation of the NLDE from first principles. In 
Sec.[3]we discuss both discrete and continuous symmetries 
common to relativistic systems, providing a clear physical 
interpretation in the present context. Finally, in Sec. [4] we 
conclude. 

2. The nonlinear Dirac equation 

2.1. Two- Component Spinor Form of the NLDE 

The second quantized Hamiltonian for a weakly interact- 
ing bosonic gas in two spatial dimensions is 



Hn 



(i) 

(2) 



The bosonic field operators -0 = ip(r, t) obey bosonic com- 
mutation relations in the Heisenberg picture. In Eq. ([1]), 
g = Anh 2 a s /m is the coupling strength for binary contact 
interactions with a s the s-wave scattering length and m the 
atomic mass. The external potential V{r) is a honeycomb 
lattice formed by standing waves of three sets of counter- 
propagating laser beams [5] . The atoms experience this po- 
tential via the AC Stark effect. We assume that the third 
spatial dimension is frozen out by a tightly confining po- 
tential which is locally harmonic, as in Ref. [5]. 

The honeycomb lattice has two sites in the lattice unit 
cell. We refer to the resulting two degenerate sublattices as 




(a) Hexagonal lattice structure 




-K 

(b) Reciprocal lattice 




(c) Nearest neighbor displace- 
ment vectors 



Fig. 1. Characterization of a honeycomb lattice. 

A and B. Expanding in terms of Bloch states in the lowest 
band belonging to A or B sites of the honeycomb lattice, as 
shown in Fig.[TJ we can break up the bosonic field operator 
into a sum over the two sublattices: 



$A = ^ a e^-^uir- r A ) 

A 

$ B = be^^-^uif - r B ) . 

B 



(3) 
(4) 

(5) 



where a and b are the time-dependent destruction operators 
at A and B sites and r A and tb are the positions of A 
and B sites, respectively. The spatial dependence is then 
encapsulated outside the operator in the exponential and 
the functions u. The summation indices indicate sums over 
A or B sites. 

Inserting Eq. ^ into Eq. ((T|), the Hamiltonian can be 
rewritten 



H= I d 2 r 



Tjj B )H (iJA + i>i 
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+ ^a+tP 1 b)Wa + ^b)(iPa+tPb)(iPa+^b)\ ■ (6) 

In the integral over Hq , imposing the restriction of nearest- 
neighbor interactions in the tight-binding, lowest band ap- 
proximation eliminates all A-A and B-B transitions except 
for on-site kinetic and potential terms; the latter can be ne- 
glected as an overall self-energy. Then only integrals involv- 
ing neighboring A-B sites remain in the sum. Similarly, in 
the interaction term only on-site terms are non-negligible, 
i.e., overlap of functions u belonging to the same site. Thus 
in the tight-binding, lowest band approximation, Eqs. (j4|- 
([5]) are substituted into Eq. ([6]) to yield: 



H= f d 2 r 



- tk -^u(x A )H ae^u(xA) 



<A,B> 

+ot e-^u(x A )H be^u(xB) 
+S* e- ii °-* B u{xB)H Q ae ii °*MXA) 
+bU- i ^ B u(x B )H be i ^ B u(xB) 

+ | J d 2 r^2a)a)aa [u(xa)} 4 
d 2 rJ2b r b r bb[u(x B )}\ 



XA 



B 

TA, XB^ 



Tb , XAB =r A -r B 



(7) 

(8) 



Here the bracketed A and B summation index signifies a 
sum over nearest-neighbor A and B sites. Isolating the in- 
tegrals by pulling out all sums and terms not dependent on 
r, Eq. ([7]) becomes 



H = 



<A,B> 



d 2 re-^ r u(x A )H e lk - r u(x B ) 



d 2 re^ k - r u(x B )H Q e^ r u(xA) 



A 



d 2 re- tk MxA)H e lk Mx A ) 



+ yyb d 2 re- tk - r u(xB)H e* k - r u(x B ) 



B 



— > a' a' aa 
2 Y 



d 2 r[u(x A )} 4 



Wbb / d 2 r[u(x B )}\ 



(9) 



Finally, we redefine the spatial integrals in Eq. © as hop- 
ping energy th and interaction energy U, respectively, as is 
standard for the Hubbard Hamiltonian 23 24] (note that 
we reserve t for time). The terms in Eq. ([9|) proportional 
to a)d and wb just count the total number of atoms in the 
system, and can been neglected as an overall constant. This 
leads to the Hamiltonian 



H = 



-th [ a 

<A.B> 



(10) 



Equation ([TO]) is the Hubbard Hamiltonian divided into 
two degenerate sublattices A and B, appropriate to the 
honeycomb optical lattice. 

In order to work towards the nonlinear Dirac equation, 
we calculate the time evolution of a and b according to the 
standard Heisenberg picture prescription. This is similar 
to the approach taken by Pitaevskii in his landmark paper 
which first obtained the NLSE, or Gross-Pitaevskii equa- 
tion [25] . The Heisenberg equation of motion is 



ihd t a k = [a k ,H] 



(11) 



The operator fife, which destroys a boson at site k on sub- 
lattice A, satisfies the bosonic commutation relation 

[6fc)Ofc/] = hk' ■ (12) 

Then the commutator with the on-site interaction terms 
reduces to 



[dk,d\d\a k a k } 



d k aldld k a k 



d\d\d k a k a k . 



(13) 



Taking the first product on the right and commuting the 
furthermost left a through according to Eq. (jT2J) , one finds: 

(14) 



a k d\d\d k d k = 2d\a k d k + a k a\a,ka k a k . 
Substituting Eq. (fT4| into Eq. (|13|) . one obtains 
[a k , a\a\a k d k ] = 2a\a k a k . 



(15) 



Substituting Eq. (fT5f into Eq. (fTTjl and the Hubbard Hamil- 
tonian Eq. (fT0|) , one finds 



ihd t d k = -tu 



b k e 



ik-{r Ak -r Bk ) 



-b k -n 2 



ik-(r Ak -r Bk 



ik-(r A 



Ud\d k d k 



(16) 



where the first three terms on the right hand side represent 
transitions from the three B-sites nearest the fc th site of the 
A sublattice and Hi and H2 are primitive cell translation 
vectors for the reciprocal lattice, as shown in Fig. [1] 

In a similar fashion as Eqs. (fTTj) - (fT6]) , we arrive at an 
expression of the same form for the B sublattice: 



ihdtb 



-t h 



a k e 



-ik-(f Ak 



.* — ik-(r A —r Bt .) 



Ub k b k b k . 



(17) 



Continuing to follow Pitaevskii's method, we next take the 
expectation value of Eqs. (fT6|) and (|17|) with respect to on- 
site coherent states. A tensor product over sites of such 
coherent states is also assumed [TIT]. A more formal, care- 
ful treatment of finite number states, rather than coherent 
states, has been worked out in the literature (see [26] and 
references therein). Either way, we obtain coupled equa- 
tions of motion for discrete, on-site, complex-valued ampli- 
tudes. For simplicity of notation we take 



a k 



(a k ) , b k = (b k ) 



(18) 
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Inserting the nearest-neighbor vectors 5%, 62, and £3 in the 
exponentials in Eqs. (j 16[) and (JTTJ) , as shown in Fig. [TJ we 
obtain 



iha k = -t h (b k e ik ' Ss + 6 fc _ ni e* Sl + b k . ri2 e lk ^) 
+Ua* k a k a k , 

ihb k = -t h (a k e- iU * + a k+ni e'^ + a k+n2 e^ 2 ) 
+Ub* k b k b k , 



(19) 
(20) 



where ri\ and n 2 in the indices label the lattice sites in the 
two directions of the primitive-cell translation vectors n\ 
and n%. 

The NLDE is derived around the linear band crossings 
between the A and B sublattices at the Brillouin zone cor- 
ners |27j , called Dirac cones in the graphene literature [3] . 
To this end, we insert particular values for the nearest- 
neighbor displacement vectors S and evaluate k at the Bril- 
louin zone corner, defined by k = K = (0,47r/3) , <5i = 
CaTs' -3). & = (aTs'D' & = (-^,0), all in the x - y 
plane. Then Eq. (fl^|) becomes 



ihh k = -t h {b k e° + b k - ni e i2jr/3 + b k _ n2 e l 
+Ua* k a k a k . 

Reducing the exponentials, 

iha k = -t h [b k + 6 fe _ ni (-1/2 - iy/3/2) 

+fe fe _„ 2 (-l/2 + iV3/2)] + Ua* k a k a k 



2tt/3 



(21) 



(22) 



In anticipation of taking the long wavelength, continuum 
limit, as is necessary to obtain the NLDE, we group terms 
appropriately in Eq. (|22|) in order to construct discrete ver- 
sions of derivatives: 

ihh k = -t h [b k + (b k - & fc _ ni )(l/2 + iy/3/2) 

-b k (l/2 + ^V3/2) + (b k ~ 6 fe _ n2 )(l/2 - iV3/2) 
-b k (l/2-iV3/2)} + Ua* k a k a k , (23) 

which reduces to 

iha k = -t h [(b k - 6 fc _ ni )(l/2 + iVa/2) 

+ (b k - 6 fe _„ 2 )(l/2 - iV3/2)] + Ua* k a k a k . (24) 

Taking the continuum limit and replacing the discrete 
quantities a k = a k (t) and b k = b k (t) by the continuous 
functions ipA = ipA(f, t) and tp B = ipB(?,t), we arrive at 



ifmpA = -t h 



d ni ip B {l/2 + iV3/2) + d n2 (l/2 - *V3/2) 



+ Ulp* A ^ A ^A ■ 



(25) 



where the partial derivatives are in the directions of the 
unit-cell vectors n± and n 2 - With a little trigonometry we 
find that the unit-cell vectors are 



n\ = cos(7r/6)e x — sm(n/6)ey = v / 3/2e x — l/2e 
H2 = cos(7r/6)e x + sin(7r/6)ej, = \/3/2e x + l/2e 



y • 



(26) 
(27) 



Up to now the "hat" symbol (circumflex accent) has been 
reserved for operators alone. However, in Eqs. (|26"]) -(|27 p we 
use this symbol to indicate a unit vector in the x and y 
directions. Derivatives with respect to the unit-cell vectors 
take the form 



d ni = nl • V = {V3/2)d x - {l/2)d y , 
d n2 = n 2 • V = (V3/2)d x + {l/2)d y . 

Substituting Eqs. Ip8| l -(|29 |) into Eq. ([25]). 
ihi> A = -t h {[{V3/2)d x - {l/2)d y ]ip B [(1/2) + i(V3/2)} 
+ {(V3/2)d x + (l/2)fy]Vn[(l/2) - »(V3/2)]} 

+ Ulp A '>pA'lpA ■ 

Further simplification of Eq. (|5D|) leads to 



(28) 
(29) 



(30) 



iHtpA — 



-(d x tpB - id y ip B ) + Utp A ipAip/ 



Similarly, for the continuum limit of b k = b k (t) 
ipB(f,t), 



ihipi 



(-d x ip A - idyil) a) + Uip B ip B ip, 



b ■ 



(31) 

IpB = 

(32) 



Equations (|3T|) - ([32)) are in fact massless Dirac equations 
with an added nonlinear term. To put this in a more fa- 
miliar form, we use Pauli matrix notation. To this end, one 
must rename the coordinate axes so that x — » y, and, in 
order to preserve the handedness of the coordinate system, 
y — > — x. We also reinsert the lattice constant a; note that a 
is unrelated to the s-wave scattering length a s briefly men- 
tioned in the definition of g following Eq. (fTJ) . Thus d x — > 
-d x . Then Eqs. pi ]) -(152 1 become 



d v and d y 



ill! • , = - th2 7r ^. [p y ^ )B _|_ id x 1p B ) + UlpA^AlpA , 



ihipi 



2 



(-dytpA + id x ipA) + Ui>* B ip B ^ B , 



(33) 
(34) 



or in matrix form, 





IpA 


— ii/jOv3 


d x — id y 




IpA 


ih 












IpB 


2 


d x + id y 




•>Pb 



-u 



^A^AlpA 
1p B 1pBlpB 



(35) 



We can write Eq. p5p more compactly in terms of Pauli 
matrices (a x , a y ) = a, 



ih 



ipA 
i>B 

Equation 



-ithCi^/3 



a- V 



4>A 
4>B 



U 



tpA^A^A 
lp* B lpBlpB 



(36) 



is the NLDE in (2+1) dimensions. 
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However, we can make one further step by expressing 
Eq. (|36p in a more covariant-looking form as follows: 



(ia d t + ic s a ■ V) 





- u 










iP* b 4>b4>b 



To this end, we seek to put the NLDE into a form con- 
sistent with the standard compact four-component spinor 
notation for the linear Dirac equation 



(37) i'fd^ = , 



(42) 



where a and V are restricted to the x — y plane. In Eq. ()37D 
thaV3 



2h 



(38) 



is an effective speed of sound of the condensate in the 
lattice [28129] (in graphene it would be replaced with the 
Fermi velocity [3]). This velocity is an effective speed of 
light for excitations of the NLDE in our QED2+1 theory. 
Experimental values of c s in BECs are on the order of cm/s, 
ten orders of magnitude slower than the speed of light in a 
vacuum. Note also that in Eq. ([37]) U has now absorbed a 
factor of l/h. 

Finally, a few additional definitions lead to a nicely com- 
pact form for the NLDE. Let 



where the 7 M are the 4x4 Dirac matrices [31] . We point out 
that this notation is not only standard but more appropri- 
ate if the quasi-particles, i.e., the long- wavelength excita- 
tions, develop a non-zero effective mass, as can be caused 
by lattice distortion [B] . 

We obtained the NLDE by evaluating the exponentials 
at the Brillouin Zone corner K + = (0,47r/3). There is 
another inequivalent corner, as shown in Fig. [2 K— = 
(0, — 47r/3), near which perturbations in momentum, i.e., 
long-wavelength quasiparticle excitations, are governed by 
a similar first-order wave equation. By considering equa- 
tions of motions derived around both Brillouin corners, we 
can obtain the four- vector notation. The coupled equations 
evaluated at /?_ are 



.4 



1 









4>A 




,B = 











1 







(39) 



With the choice of metric which raises and lowers space- 
time indices restricted to (2+1) dimensions, 



ihipA = th< ^ (d x + id y )ip B + Uip* A ip A ip A , 

ihipB = —^Y^-{~d x + id y )i/j A + U^* B i/j B ipB ■ 
Following the same steps as before we obtain 



ih 



1 

0-10 
00-1 



IpB 
4>A 



-a- V 



tpB 
IpA 



u 



l^B^BlpB 
fpA^AlpA 



(43) 
(44) 

(45) 



(40) 



We combine Eqs. (|45[) and (|36[) into one equation involving 
a single 4-component object and attach ± subscripts to 
the wave functions '5 to specify the corner of the Brillouin 
Zone. The resulting expression is 



Eq. ([37]) becomes 




*+ 






id t 


+ ic s 


{ia^df, - UtpAijjA ~ U^BifiB)^ = , 


(41) 







S-V 






-a- V 







where the standard Einstein summation rule is in effect, 
/i € {0,1,2} in keeping with (2+1) dimensions, and the 
units are chosen such that c, = 1. 



2.2. Maximally Compact Form of the NLDE 



-U 



= 0. 



(46) 



where the nonlinear terms are grouped into the 2-vectors 
N + and iV_ defined by 



In Sec. 12.11 we developed ip, a two-dimensional complex 
object which brings to mind one member of a pair of Weyl- 
spinors in the (1/2,1/2) chiral representation of the Dirac 
algebra used to describe massless neutrinos in the stan- 
dard QED3 + i theory. Such a treatment is appropriate for 
any neutral Dirac fermion viewed in the extreme relativis- 
tic frame. In order to make the connection clear we must 
find the second member of the pair of Weyl-spinors and ver- 
ify that the mapping is true. A thorough treatment of the 
mapping of QED3 + i into the QED2+1 theory of graphene 
is contained in [30 and references therein. We will continue 
to restrict ourselves to (2+1) dimensions in the following. 



7V+ = 



iV_ = 



(ip B 1p B 1pB)- 
(lp A 1p A ^A)- 



(iPa^a^Pa)+ 
(V>bV>bV , b)+ 
and the four-spinors are given by 

lpA+ 

l\>B- 
IpA- 

^A+ ^*B+ IpB— *1>A- 



(47) 



\Jr* vj/* 



(48) 



(49) 
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Here again the ± subscripts refer to the specific corner of 
the Brillouin zone. 

We can reduce Eq. (|46|) to a more compact form by in- 
troducing the following 4x4 matrices: 





1 




a x 






,Si = 






1 




-(J 




A 







A+ = 




,A- = 











A 




B 







B+ = 




,B- = 











B 





-<7„ 



(50) 
(51) 
(52) 



The boldface notation A and B denote the 2x2 matrices 
defined in Eq. (|39|) . Also, the boldface entries 1 and in 
the matrices in Eq. (|47l) refer to the 2x2 unit matrix and 
zero matrix, respectively [15] . Then Eq. (|46|) becomes 



- U^A-VA- - UV f B-V B-)V = . 



(53) 



We substitute into Eq. (f5"3"]) the Dirac matrices in the Chiral 
representation: 



1 




-a x 




-<Ty 








,7 2 = 




1 




a x 




_a y _ 



and multiply on the left by 7 . Note that /x € {0,1,2} in 
keeping with (2+1) dimensions. Finally, for the nonlinear 
term we introduce 



N : 



E 



(55) 



Implementing Dirac matrix notation as described, the 
NLDE of Eq. {53]) becomes 



(ry^ + 7°iV)* = 0. 



(56) 



Equation (fB"6"|) is the most compact form of the massless 
NLDE in 4-component spinor notation. 

3. Symmetries and Constraints 

The NLDE as expressed by Eq. ([37]) or Eq. (jUJ) looks 
like the Dirac equation for a two-spinor, as stated in the 
graphene problem 3 , with the addition of two nonlinear 
on-site interaction terms, one for each A and B sublattice; 
similarly, the more compact form developed in Sec. l2.2l also 
appears to be a Dirac equation for a four-component spinor, 
with an additional term. However, one should not be too 
hasty in assigning characteristics based on appearances. We 
therefore make a careful and thorough exploration of the 
symmetries and other important mathematical properties 
of the NLDE. 



In what follows we follow a similar route as in Ref. [TS] . 
First we check the linear version with the delta interactions 
turned off to ensure that Eq. (|4T]) is indeed the massless 
Dirac equation in (2+1) dimensions with all the necessary 
symmetries. For each symmetry we then check the nonlin- 
ear interaction terms to determine whether they preserve 
or break the symmetry. 

3.1. Locality 

We do not necessarily require the evolution of the wave- 
function as described by an NLDE to be governed by a lo- 
cal theory. A local theory is one in which the terms in the 
linear equations of motion involve only factors of the wave- 
function and its derivatives evaluated at the same space- 
time coordinate. Nonlocality arises in low energy limits of 
some quantum field theories. However, the nonlinearity in 
our NLDE is manifestly local. Thus our NLDE is closer to 
the standard Dirac equation on the classical level (no quan- 
tum effects), modified by the on-site interaction term. In 
Sec.|4]we discuss the possibility of non-local nonlinearities, 
including for graphene. 

3.2. Poincare Symmetry 



A Poincare transformation takes the spatiotemporal 
(54) point defined by the the 4-vector r v into the point r' 



according to 



(57) 



where A^ is the coordinate matrix representation of the 
Lorentz group and is a space-time translation. The wave 
function ^ transforms as 

= M(A)*(r) 

where the matrices M(A) form a representation of the sub- 
group of the Lorentz group consisting of spatial rotations 
and boosts. Boosts can be thought of as rotations in imag- 
inary time by imaginary angles mixing space and time co- 
ordinates. We restrict these transformations to (2+1) di- 
mensions. 

The proof of Lorentz covariance of the standard massless 
Dirac equation, Eq. (j42]) . is arrived at with the aid of the 
transformations for the wave function and partial deriva- 
tives: 



*(r) = M- 1 (A)*'(r') 

d x \?)y . 



(58) 
(59) 



This yields the conditions for the Dirac matrices: 

7j = AjiMfiM- 1 . (60) 

The standard form of the Dirac matrices obtained this way 
can be found in the literature |31] and are identical to the 
results of our theory. 
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Thus imposing Lorentz covariance on the NLDE as 
expressed in Eq. (|56[) requires \& to transform under 
the irreducible representation of a subgroup of SL(2,C), 
the 2x2 complex matrices of unit determinant. The 
four-dimensional representation of SL(2,C), Z^ 1 / 2 ' 1 / 2 ' , is 
formed by taking the direct product of the two-dimensional 
representations D^/^and D^ 1 / 2 ); 



£,(1/2,1/2) = £,(1/2,0) (g, £,(0,1/2) 



(61) 



This subgroup of SL(2,C) is isomorphic to a subgroup of 
the Lorentz group, the one obtained by restricting Lorentz 
transformations to the plane of the honeycomb lattice. So, 
the upper two components of \& transform as a spinor under 
£)(i/2,o)^ foe lower two as a spinor under I)( ' 1 / 2 ) ; and ^ 
itself as a four-component spinor or bispinor |32j . 

The next task is to examine the behavior of 'J, as defined 
in Eq. (|49p and governed by Eq. (|56p . under rotations in 
the x — y plane. In order to obtain Poincare covariance we 
must show that it is the same as that of a four-component 
spinor in the standard Dirac3 + i theory restricted to the 2D 
plane. The honeycomb lattice is invariant under rotations 
by ±27r/3 but the four components of \1/ are also defined 
by the particular corner of the Brillouin Zone. Since we're 
considering a discrete lattice it is natural to discuss discrete 
rotations which realign lattice points and in the continuum 
limit map the continuous rotations of QED2+1 onto our 
theory. Rotations by ±7r/3 exchange A and B sites, and 
take the theory to that of the opposite K point: K+ does 
not go to K- but the result after calculating the relative 
phase exponentials gives back the same theory. To see this 
we chose a different primitive unit cell, the one obtained 
by a rotation of 27r/3 about the n\,ni origin in Fig[TJ This 
is because the direction of K\ (defined as K + rotated by 
7r/3) differs from that of by 27r/3. 

Thus under this discrete rotation 



b k - 


* bk-2nx , 


(62) 




* bk—ni —TI2 5 


(63) 


bk — U2 


* bk—ni ! 


(64) 


au - 


* dk—ni j 


(65) 


dfe+ni — 


* Q j k — 2n\ +712 1 


(66) 


<2fc+n 2 ~~ 


* a/c-2ni ■ 


(67) 


Also we observe that 




K'+ • S 3 


= • Si 


(68) 


K'+ ■ A 


= K--S 2 


(69) 


K'+ • 5 2 


= K--5 3 


(70) 



Substituting Eqs. j62|)-([70]) into Eqs. (Tl9])-(|20|) we obtain 

iha k - ni = -th(bk-2m e lK ~' 1 + &fc_ ni _„ 2 e lK -' S2 

+b k - ni e lK - •*> ) + Ua* k _ ni a k - ni a k - ni , (71) 



Now we redefine the index k: in Eq. (fTTj) . k 
in Eq. O, k k + 2n\. Then 



) + Ub* k _ 2ni bk-2mb k -2nf (72) 
k + ni, while 



iha k = -th(bk-m e lK ' 5l + 6 fe _„, 
+b k e iS -- Ss ) + Ua* k a k a k 

ihb k = -t h (a k+n2 e ~ lK -' Sl + a k+ 
+a k e^ R - S3 ) + Ublb k b k 



lK--8 2 



(73) 
(74) 



To summarize, a rotation by 7r/3 which takes 
is identical to the unrotated theory but with K + — > 
Note that the redefinition of the index k is different for the 
two equations because old and new primitive cells are re- 
lated by a rotation. Rotating by ±27r/3 exchanges A and B 
sites once more and returns ^ to its original configuration. 
Since ^ has four components, the effect of making one full 
rotation of 2tt is that the components acquire a net phase 
so that ^ — > — This Berry phase [2] endows \& with the 
characteristic double- valuedness of a genuine 4-spinor. 

However, we must be cautious when discussing chirality 
and helicity, since we treat (2+1) dimensions and can use 
only the first two Pauli matrices. As in the (3+1) theory, one 
can define a pseudo-chirality operator in (2+1) dimensions, 
7 5 , as the product of the other four 7 matrices. In the Weyl, 
or Chiral representation we have 



7 5 = i7°7 1 7 2 7 3 



1 
1 



(75) 



where again the boldface indicates a 2 x 2 submatrix. This 
is the natural representation for *5> in that the NLDE of 
Eq. (|56p maps into this representation in a natural way: 
states of well defined chirality correspond to the upper and 
lower 2-spinors. Thus the upper and lower spinors 



(76) 



v+ 



















are eigenfunctions of 7° , 
7 5 *+ = , 7 5 *_ = 



(77) 



The question of Poincare covariance of the nonlinear 
Dirac equation remains. First, we check coordinate trans- 
lations. The wave function is required to transform as 



*'(r^ + d„) = *(r p ). 



(78) 



ihb k -2 ni — —th(o>k-m 



e + ak-2 ni +n 2 e 



Thus the nonlinear terms which have two factors of are 
invariant under translations. Under spatial rotations we ob- 
serve that the interaction terms in the NLDE remain un- 
changed within the context of the full theory. We include 
both K + and K- points in the full theory. As for the case 
of boosts, we note that for the linear equation the compo- 
nents of the wave function transform in accordance with 
the transformation of the space-time coordinates in their 
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arguments. This is exactly canceled by the reciprocal trans- 
formations of the partial derivatives. This is not the case 
for the nonlinear terms. One has some matrix product with 
two factors of the wave function. Thus the nonlinear Dirac 
equation is not invariant under Lorentz boosts. 

For any truly fundamental theory of nature it is de- 
manded that the governing equations be invariant with re- 
spect to the Poincare group. Sometimes this requirement 
is loosened as in the quantization of gauge theories, when 
fixing the gauge causes relativistic invariance to be non- 
manifest. Yet the theory itself certainly remains invariant. 
Other times the breaking of Poincare covariance implies 
that we are dealing with an effective theory in which deeper 
physical processes are at work. Our case is an example of 
the latter. By including on-site interactions we break rela- 
tivistic covariance of the linear equation of motion by in- 
troducing self-interactions which can be viewed as evidence 
of "deeper physics" in our two-dimensional universe. 

3.3. Hermiticity 

We require the Hamiltonian be Hermitian in order to 
guarantee that physically measurable quantities (observ- 
ables) are real. Thus each term must be independently Her- 
mitian. In particular, we must show that — N. The 
proof follows: 

Art = (u&A+VA+y 
= UA\.&A\.(&)i 
= UA + &A + V 

= N (79) 

where the last two steps work because A + is real and sym- 
metric. The nonlinear terms, and therefore the NLDE, are 
indeed Hermitian. 

3.4. Current Conservation 

Current conservation is expected in ordinary QED for a 
closed system, i.e., an isolated volume of space. However, 
most theories of quantum gravity do introduce Lorentz vio- 
lating terms which bring in charge non-conservation effects. 
The conserved current for the linear Dirac equation is 



y 



(80) 



where iff = ^f'^ 

We check that current is also conserved for the NLDE. 
The 3-divergence of the current in (2+1) spatial dimensions 
is 



where £ {0,1,2} and we have used the chain rule and the 
anti-commuting properties of the Dirac matrices. Taking 
the adjoint of the NLDE in four-component form, Eq. ([56]), 
yields 



which implies 



(82) 
(83) 



Then the (2+l)-divergence of the current is 

d^f = -i^A^V* + * t 7°(i7°A r *) (84) 



-i* t (iV t - AO* = 



(85) 



Thus current is conserved. This is in fact the statement 
that hermiticity implies current conservation. 

3.5. Chiral Current 

The conserved chiral current for the linear Dirac equation 

is 

% = * 7 M 75* ■ (86) 
Then the (2+l)-divergence of the chiral current is 

« = ^(*7*75*)- (87) 
Following a similar route as in Sec. 13. 4i we find 

M = -^(^75 + 75 AO* (88) 

where we have used the hermiticity of N. This means that 
in order for chiral current to be conserved we must have 



{^,75} = 



(89) 



where the curly braces signify the anti-commutator. The 
anti-commutator is 

{N+, 75} = U^A+^A+'ys + j5U^A+^fA+ . (90) 
Writing this out in explicit matrix form we find that 



{N +>l5 } = +2N, 
{N_, l5 } = -2N. 



(91) 
(92) 



(81) 



Therefore chiral current is not conserved. This fact is remi- 
niscent of the anomalous non-conservation of chiral current 
in the case of some field theories upon quantization, which 
is mediated by instantons. It is interesting to consider that 
our nonlinear terms might be treated as quantum-induced 
nonlinearities. 

3.6. Universality 

Universality refers to the invariance of a theory under 
rescaling of the solution. For the case of the linear Dirac 
equation, * — > A* leaves the theory unchanged. Because 
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the nonlinear term in the NLDE contains one factor of ^ 
and one factor of its adjoint, it scales as A 2 , thus breaking 
the invariance of the nonlinear theory. For a treatment of 
NLDEs which are universal in this sense, but not relevant 
to the present solid state system, see Ref. [15] . 

3.7. Discrete Symmetries 
3.7.1. Parity 

In the case of the standard massless Dirac equation, in- 
variance under the parity operation P requires \& to trans- 
form as 



= = 7°*, 



(93) 



where, in contrast to the Chiral representation presented 
in Eq. ([54]) . we take 7 in this section to be in the Dirac 
representation, 



1 
1 



(94) 



We proceed to determine if the NLDE remains invariant 
under the transformation of Eq. (193|) . The parity operator 
acting on the honeycomb lattice inverts both coordinate 
axes, V — > —V, and thus exchanges A and B sites while 
also exchanging if-point indices. The transformed linear 
equations (U = 0) are 







-a- V 


7° 








+ ic s 






= 


.*+. 




& • V 









(95) 



Interchanging upper and lower spinors we obtain the equiv- 
alent form 









i 




+ ic s 









a ■ V 

— 5 • V 



(96) 



Thus the linear equations are invariant under parity. 

Next we check the nonlinear term in the NLDE. With the 
transformed nonlinear term added to Eq. (|95p we obtain 









% 




+ic s 









-CT-V 

a- V 



-U 



= 0. 



(97) 



which after similar steps as above gives us back the NLDE. 
Therefore parity is a symmetry of the nonlinear Dirac equa- 
tion. We note that parity is conserved in standard QED 
and also in the theory of the strong interactions, but vio- 
lated by weak interactions. Therefore, in this aspect, the 
NLDE is consistent with QED. 



3.7.2. Charge Conjugation 

In order to maintain symmetry under charge conjugation 
one requires the nonlinear term to transform as 



( 7 °A0'= Cj°N*C- 



(98) 



where in the Chiral representation the charge conjugation 
operator is 

C = l l 

and the wavefunction transforms in the standard way: 



7 7 



(99) 



We determine if the NLDE maintains this symmetry. Let 
F be one of the nonlinear terms. Then 

( 7 °JV)' = ( 7 1 7 * tT ) t ^+(7 1 7°* tT M+ 

= f T 7 ° 7 1 A + 7 1 7°**A + (100) 

Also, we have 

CfN*^ 1 =7 1 7°(* t ^+*^+)*7 1 

= 7 1 7°* T ^+**-4+7 1 (101) 
Thus the NLDE breaks charge conjugation symmetry. 



3.7.3. Time Reversal 

The usual time reversal, t 
form as 



* -> = 6* = i 7 1 7 3 *, 
where O is the time-reversal operator and 

-<7,, 



t, requires that \& trans- 
(102) 



i7 X 7 3 = 







(103) 



In our theory the intrinsic effect of time reversal is to change 
the direction of momentum so that K points are switched 
without exchanging A and B indices. Combining these ef- 
fects, we determine whether or not the linear part of the 
NLDE, i.e., the linear Dirac equation, remains invariant: 



■ <5 • 1 3 



+ic s 

Then 

. d 
- l di 



a- V 
















+ ic s 











Z7 1 7 3 






.*+. 


a- V 








a- V 



= 



(104) 



= 



(105) 



The appearance of the negative sign in front of the time 
derivative indicates that these are the equations satisfied 



by the negative energy solutions, or holes. Thus time rever- 
sal takes the equations describing electrons into holes and 
those for holes into electrons but keeps the overall theory 
invariant. 

We proceed to consider the nonlinear term. With the 
transformed nonlinear term we obtain 









+ ic s 






N + 






= 0. 



A_ 

In Eq. (|106p the interaction term has acquired a factor of i. 
Therefore the NLDE is not invariant under time reversal. 
Thus, as one also observes in the Standard Model of particle 
physics, CP and T are not conserved independently but 
CPT is conserved. 

4. Discussion and Conclusions 



be looked for in quantum electrodynamics under the proper 
circumstances [21] , 

All generalizations of the scalar nonlinear Schrodinger 
equation relevant to BECs are candidates for generalized 
nonlinear Dirac equations. For example, pseudo-spin struc- 
ture leads to a vector NLSE. One can therefore anticipate 
a vector NLDE as well. 

We have thoroughly explored both continuous and dis- 
crete symmetries of the NLDE. In particular, we showed 
that pseudo-chiral current is not conserved, the NLDE is 
not covariant under Lorentz boosts, and it breaks charge 
conjugation as well as time reversal symmetry. On the other 
hand, the NLDE is hermitian, local, conserves current, and 
is symmetric under parity. The NLDE maintains CPT sym- 
metry, just as one finds in the Standard Model. 

In a future work we will treat soliton and vortex solutions 
of the NLDE, as a first step towards a complete classifica- 
tion of nonlinear rclativistic phenomena in BECs. 

We acknowledge useful discussions with Alex Flournoy 
and Mark Lusk. This work was supported by the National 
Science Foundation under Grant PHY-0547845 as part of 
the NSF CAREER program. 



The nonlinear Dirac equation we have presented intro- 
duces a completely new class of nonlinear phenomena in 
Bose-Einstein condensates. Although our work is related 
to graphene, in that the BEC is taken to be trapped on a 
honeycomb lattice, we have switched bosons for fcrmions. 
The form of the nonlinearity is then a natural physical re- 
sult of binary interactions between bosons. In fact, it is a 
spinor generalization of the kind of nonlinearity one finds 
in the nonlinear Schrodinger equation, i.e., proportional to 
the local condensate density. The same equation will occur 
for light subject to a Kerr nonlinearity when propagating 
through a photonic crystal with a honeycomb lattice struc- 
ture [33] , 

We showed that the NLDE breaks Poincare covariance, 
and therefore the principle of relativity. This suggests that 
small nonlinearities of this form could be looked for by such 
symmetry breaking in a variety of systems where Dirac or 
Dirac-like equations apply. For instance, even for fermions 
there is a small mean field effect. We could just as easily 
have considered ultra-cold fermions on an optical lattice. 
This would appear at first sight to be the exact analog of 
graphene; however, our work shows that Poincare covari- 
ance will be broken by mean field effects, even if on a small 
level. Indeed, for graphene one should expect similar ef- 
fects due to Coulomb interactions. The latter nonlinear- 
ity can be expected to be non-local due to the power law 
behavior of 1/r for the Coulomb potential, just as dipole- 
dipole interactions between ultra-cold atoms even without 
an optical lattice lead to a non-local nonlinear Schrodinger 
equation in a 3D continuum (in 2D dipole-dipole interac- 
tions, which have a power law of 1/r 3 , are local). Thus, in 
graphene, we make the conjecture that there is a non-local 
nonlinear correction term to the massless Dirac equation 
which breaks Poincare covariance. Similar corrections can 
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